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Abstract 

Solving the 4-d Einstein equations as evolution in time requires solving equations of two types: 
the four elliptic initial data (constraint) equations, followed by the six second order evolution 



o 

I equations. Analytically the constraint equations remain solved under the action of the evolution, 



and one approach is to simply monitor them (unconstrained evolution). 

The problem of the 3-d computational simulation of even a single isolated vacuum black hole 
has proven to be remarkably difficult. Recently, we have become aware of two publications that 

(N 

J> ■ describe very long term evolution, at least for single isolated black holes. An essential feature in 

. each of these results is constraint subtraction. Additionally, each of these approaches is based on 

I what we call "modern," hyperbolic formulations of the Einstein equations. It is generally assumed, 

' based on computational experience, that the use of such modern formulations is essential for long- 

O . term black hole stability. We report here on comparable lifetime results based on the much simpler 

. ("traditional") g - K formulation. 

With specific subtraction of constraints, with a simple analytic gauge, with very simple boundary 
^ ■ conditions, and for moderately large domains with moderately fine resolution, we find computa- 

H : 

^ tional evolutions of isolated nonspinning black holes for times exceeding lOOOGM/c^. 

We have also carried out a series of constrained 3-d evolutions of single isolated black holes. We 
find that constraint solution can produce substantially stabilized long-term single hole evolutions. 
However, we have found that for large domains, neither constraint-subtracted nor constrained g 
- K evolutions carried out in Cartesian coordinates admit arbitrarily long-lived simulations. The 
failure appears to arise from features at the inner excision boundary; the behavior does generally 
improve with resolution. 

PACS numbers: 
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I. INTRODUCTION 



Binary black hole systems are expected to be the strongest possible astrophysical gravi- 
tational wave sources. In the final moments of stellar mass black hole inspiral, the radiation 
will be detectable in the current (LIGO-class) detectors. If the total binary mass is of the 
order of lOM©, the moment of final plunge to coalescence will emit a signal detectable by 
the current generation of detectors from very distant (Gpc) sources. The merger of super- 
massive black holes in the center of galaxies will be the dominant signal in the spaceborne 
LISA detector, and detectable out to large redshift. 

Simulation of these mergers will play an important part in the prediction, detection, and 
the analysis of their gravitational signals in gravitational wave detectors. To do so requires 
a correct formalism which does not generate spurious singularities during the attempted 
simulation. Recent important work^j^l has been done in extending the computational life- 
time of single isolated black hole simulations. We report here on such an extension which 
demonstrates that constraint subtraction by itself is adequate to produce very long-lived 
simulations. We demonstrate this even for a very simple ( "traditional" ) g - K formulation, 
with specific subtraction of constraints (which are analytically zero) for single isolated non- 
spinning black holes, with a simple analytic gauge (lapse and shift not "densitized" ; see 
IIVAI below), with very simple boundary conditions, and for moderately large domains with 
moderately fine resolution. 

We have also found that constrained 3-d evolution with densitized lapse can produce 
substantially stabilized long-term single holes, even for subtractions that differ from the 
values that we have found to be optimal in the unconstrained case, and we give some 
preliminary constrained evolution results. 

However, in all cases we find that attempting to carry out the simulations on very large 
domains (±20M, or larger) still yields eventually unstable simulations, by any of the methods 
reported here. 

II. 2, + I FORMULATION OF EINSTEIN EQUATIONS 

:e a Cauchy formulation (3+1) of the ADM type, after Arnowitt, Deser, and 
. In such a method the 3-metric gij and its momentum Kij are specified at one 
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initial time on a spacelike hypersurface, and evolved into the future. The ADM metric is 

ds^ = -{a^- + 2(3i dt dx' + gtj dx' dx^ (1) 

where a is the lapse function and /3* is the shift 3-vector; these gauge functions encode the 
coordinatization.f 

The Einstein field equations contain both hyperbolic evolution equations and elliptic 
constraint equations. The constraint equations for vacuum in the ADM decomposition are: 

H =^[R-K,,K'^ + K^] = 0, (2) 

W = Vj (K'^ - g'^K) = 0. (3) 

Here R is the 3-d Ricci scalar constructed from the 3-metric, and Vj is the 3-d covariant 
derivative compatible with gij. Initial data must satisfy these constraint equations; one may 
not freely specify all components of gij and Kij. 

In this paper we are concerned only with single isolated black holes. From this point 
of view the problem is not solving the initial value equations, since the data are known 
analytically. Instead, the question is one of the stability of the solution as these data are 
evolved computationally. The evolution equations from the Einstein system are 

gij = -2aKij + V, A + "^iPj (4) 

and 

K., = - V, V,« + «(% - 2K,kK) + KK,,) + P'VkK,, + K,kV,(3' + Kj^Vi^" (5) 

where a dot ( ' ) denotes the partial derivative with respect to time, and Rij is the 3-d Ricci 
tensor. 

We call this form of the Einstein equations of ADM type, referring to the fundamental 
development 0]; this specific form is called the g - K form. Here, Eq. Q-Eq. (jH)), the 
constraint equations, are the vacuum Einstein equations ^Goo = and ^Goi = respectively. 

f] Latin indices run 1, 2, 3 and are lowered and raised by gij and its 3-d inverse g^^ . 
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Eq. (@I)-Eq. (0), the evolution equations, are a first order form of the vacuum Einstein 
equations ^Rij = 0. 

The true ADM form writes the evolution equations as "^Gij = 0, the space components 
of the 4-d Einstein tensor, rather than the Ricci tensor. Frittelli and ReulaQ] have shown 
that with certain (rather strong) assumptions, there is stable maintenance of the constraints 
under unconstrained evolution for Eq. (@I)-Eq. but only neutral stability for ^Gjj = 0. 

III. DATA FORM 

In this paper we consider only single isolated black holes, so the data setting problem 
is already solved; we use Kerr-Schild data, which describes a single isolated spinning or 
nonspinning black hole hole. 

The Kerr-Schild form of a black hole solution describes the spacetime of a single black 
hole with mass, m, and specific angular momentum, a = j/m, in a coordinate system that 
is well behaved at the horizon: 



where 77^1, is the metric of fiat space, Hks is a scalar function of x^, and is an (ingoing) 
null vector, null with respect to both the background and the full metric. 



Comparing the Kerr-Schild metric with the ADM decomposition Eq. (^, we find that 
the t = constant 3-space metric is: Qij = 6ij + 2HKsh^- Further, the ADM gauge variables 
are 



ds^ = T]^^ dx^ dx'^ + 2Hks{x'')1^1^ dx" da:^ 



(6) 



0. 



(7) 



(8) 



and 



1 



(9) 




(10) 
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Each term on the right hand side of this equation is known analytically; in particular, 
for a black hole at rest, cjij = 0. 

The general non-moving black hole metric in Kerr-Schild form (written in Kerr's original 
rectangular coordinates) has 

-^^•5 = — 2 2^' ^^^^ 

+ cos"' a 

and 

_ / rx + ay ry - ax z\ 
V'' r^ + a^^ r^ + a^' rj' 

where r, 6 (and 0) are auxiliary spheroidal coordinates, z = r{x,y, z) cos6, and is the 
axial angle. r{x,y,z) is obtained from the relation, 

^ + 4 = 1. (13) 

giving 



1. o 1 



r 



2 / „2 „2 



2" ' V 4 

with 



(p2-a2) + W-(p2_a2)2 + a2z^ (14) 



p= ^x2+?/2 + ;22_ (15) 

In the nonspinning case, one has li = Xi/r, so that a = ^^^^^ ^ , and /5j = 2mxi/r. 



IV. CONSTRAINT SUBTRACTION 

The difference between Eq. (HJ-Eq. (0), and ^Gjj = 0, is a specific subtraction of the 
constraint equations. This has led to the consideration by a number of groups, of constraint 
subtraction with coefficients chosen by numerical search, or by analytical estimate (perhaps 
combined with numerical search) to improve the long-term stability of the unconstrained 
evolution. We have carried out such a numerical search, and we use the following constraint 
subtraction: 



- aH{0A64: gij + 0.36 Kij ) (16) 

on the right hand side of the Kij equation (Eq. (0)). We have found that this subtraction 
substantially improves the unconstrained evolution of nonspinning single-hole data. For 
these evolutions we choose fixed (Dirichlet) outer boundary conditions set equal to the 
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analytical value. In every case the simulation excises the interior of the black hole. One- 
sided differencing is used near the inner (mask) boundary, so no boundary condition is needed 
there, consistent with the property of the horizon as a causal boundary. The mask is specified 
at a radius of 0.75 M. The solution employs the SUNDIALS package for time integration 
0,3]. The spatial discretization is fourth order, and the time evolution is typically fourth 
order (variable-order according to the relative and absolute error tolerances). This approach 
successfully stabilizes the evolutions for domains of ±10M, as shown in Figure ^ For larger 
domains, however, the system is increasingly shorter lived, as shown in Figure |21 As a simple 
measure of the quality of solution, we present the either the I2 norm or rms norm of the 
Hamiltonian constraint constructed via a straightforward fourth order differencing scheme 
from the code results. 

Although it is difficult to compare subtraction techniques across different formal represen- 
tations of the Einstein equations, we do find typically much smaller coefficients of subtraction 
(of order 0.5) than found for different formulations, e.g. ^ with a constraint subtraction of 
order —12. The BSSN formulation of ^ is also substantially different from ours (there is a 
subtraction from the Qij equation, for instance, which we do not have, and the subtraction 
from the Kij equation is different from ours), though the coefficient of subtraction from Kij 
for this approach is small, comparable to ours. 



A. Densitized Lapse 

There is extensive evidence in the literature that a densitized lapse improves the 
hyperbolicityy] of the (at best) weakly hyperbolic ADM form of the Einstein equations. 
We implement densitized lapse for single black hole simulations by writing 

a = aanalytic(fi'/fi'analytic)^, (17) 

where aanaiytic is the explicit lapse as a function of coordinates given by Eq. S'anaiytic is the 
analytic Kerr-Schild 3— metric determinant as a function of coordinates, (gfanaiytic = l + 2i?/f , 
Ql), g is the computational 3— metric determinant, and p is an adjustable positive constant 
usually taken to be | or |. In fact we find that densitized lapse does not enhance constraint- 
subtracted lifetime, though it does contribute substantially to longer lifetime in constrained 
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evolutions described below. 



V. CONSTRAINED EVOLUTION 

The evolution of Kerr-Schild data must continue to satisfy the constraint equations, 
Eqs. (0)-© as we evolve away from the initial data. However, if we choose nonoptimal 
(e.g. zero) constraint subtraction for even nonspinning black holes, the evolution leads to 
an eventual violation of the constraints. Hence we have investigated constrained evolution, 
solving the constraint equations as part of the time update of the evolution equations. 

The postevaluated tracking of constraint errors [residual postevaluation] shown in Figures 
^Oused the direct discretization of the constraint equations Eq. Q and Eq. 0. For con- 
strained evolution we need instead to implement an accurate, efficient method of constraint 
solution. We adopt the conformal transverse-traceless method of York and collaborators llOl l- 
|14| which consists of a conformal decomposition with a scalar that adjusts the metric, 
and a vector potential that adjusts the longitudinal components of the extrinsic curva- 
ture. The constraint equations are then solved for these new quantities 0, such that the 
complete solution fully satisfies the constraints. 

Applying this approach to constrained evolution, the metric and traceless extrinsic cur- 
vature in the middle of a timestep (after an explicit integration forward in time) are taken 
as conformal trial functions cjij and A^^ . 

The physical metric at the end of the full timestep (i.e. after the constraint equation 
solve), Qij, and the trace-free part of the extrinsic curvature at the end of the full timestep, 
Aij, are related to the background fields through a conformal factor: 



Here </> is the conformal factor, and (Iw) will be used to cancel any possible longitudinal 
contribution, w"^ is a vector potential, and 




9ij = <f> 9ij, 



(18) 
(19) 



{Iw) 




(20) 



The trace K is not corrected: 



K 



K. 



(21) 
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Writing the Hamiltonian and momentum constraint equations in terms of the quantities in 
Eqs. ffTS )) - (PT|) . we obtain four coupled eUiptic equations for the fields and [l^ : 

2 ~. 



(1/8) [R<P + ^K'<P' - 



(22) 
(23) 



These equations are solved to complete each time-update step. The resulting solved Qij 
and Kij are taken as the data for the next time-update. Notice that these equations require 
no specific gauge choice. A similar approach also can be applied to other formulations which 
generally have a larger number of constraints. 



A. Elliptic Equation Boundary Conditions 

A solution of the elliptic constraint equations requires that boundary data be specified 
on both the outer boundary and on the surfaces of any masked regions. For the elliptic 
solution here we can choose simple conditions, 0=1 and = 0, on the masked region 
surrounding the singularity. Because we solve the problem on a finite domain, we also must 
provide an outer boundary condition for (p and w^. For this demonstration of the technique, 
we choose the same conditions at the outer boundary of the domain: = 1 and = 0. In 
long term evolution we expect the evolved solution to converge (as the solution is refined) 
to a solution of the constraints, so a global solution = 1 and ty* = is expected in this 
analytic limit. For achievable resolutions, however, the quantities and deviate from 
this prediction. 



VI. CONSTRAINED EVOLUTION RESULTS 

The elliptic constraint equations are solved either by a PETSc |l3, Q, Q GMRES 
solver or KINSOL GMRES solver; the spatial differencing is fourth order. We present 
below (Figures El- EI) the results of preliminary constrained evolution of nonspinning black 
holes. Figure 121 shows the rms norm of Hamiltonian and momentum constraints for these 
simulations. Compared to the relatively short term crash of the unconstrained evolution, 
the constrained evolution clearly does stabilize single black hole evolutions in small domains, 
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regardless of the precise subtraction. To fully understand the content of Figure El consider 
that, even if the residual limit in the solution of the constraint equations, Eq. (j^^ - Eq. ((221); 
is set to extremely small values (it can be set very near to machine precision, meaning that 
the discretized matrix form of equations Eq. (j22|) - Eq. (j23p can be solved to fractional 
errors of order 10^^^), the post-evaluations of the constraint residuals will typically show 
the expected zero only to the internal discretization accuracy (here fourth order). This 
is why the constrained solution shown in Figure El shows a finite (but convergent) level of 
Hamiltonian constraint violation. 

Figures mini show a constrained g — K evolution with densitized lapse, at three resolutions 
(M/5,M/7.5,M/10). The M/5 evolution became unstable before t = lOOM. The M/7.5 
run began showing large residuals at t ^ 350M. The M/10 run shows better behavior than 
M/7.5, at least initially It shows a similar late time instability which tracks (at smaller 
error) the behavior of the M/7.5 case, but around 150M ceases to be convergent; see Figure 
El Figure El shows the 2d 2 = behavior of the residual component Gxx at t = lOOM. The 
"red-blue" pattern of the features near the excision mask indicate that most of the error 
develops there. The residual becomes more asymmetrical at later times. 



VII. CORRECTNESS OF CONSTRAINED EVOLUTION 



The constraint maintenance approach uses what has been called in magnetohydro dynam- 
ics, a projection methodj2^. This method for constrained evolution raises questions about 
the meaning of the solutions obtained. This is sometimes put bluntly: "Accepting that the 
method finds solutions of the full Einstein system, how do we know that the found solution 
is the right one?" By this is meant that the constraint solution step may somehow move the 
solution back to an "erroneous" point on the space of constraint solutions. For instance, 
it might be possible that although the evolution substep and the constraint substep are 
individually convergent computational processes, the result of combining them is in some 
manner not convergent. (In the much simpler MHD case there are analytical proofs that 
projection minimizes the resultant error in the magnetic field, in a convergent way. [2^) 

There are several parts to the response. (It will be clear that we do not pretend to 
a rigorous analytical proof.) To begin with, we have constructed completely independent 



"residual evaluators" for the full Einstein system 



22| |. These evaluate the Einstein tensor. 
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working just from the metric produced by the computational solution. They are completely 
different from the way the equations are expressed in the constrained evolution code. As we 
show in Figures HEl the resulting residual is in every case initially small (order of truncation 
error) and convergent. Thus we have achieved a computational solution to the Einstein 
system. Note that our full Einstein equation residual evaluator checks both the constraints 
(Einstein equations at one time), and the evolution equations connecting different time steps. 

The residual evaluators are written to return fourth order accurate results. Since they 
have now been verified and show convergence of the solution, we appeal to the assumption 
that the Einstein system is not singular at the solution manifold. Thus we expect that the 
computational result converges to the analytical solution of the Einstein equations. We 
have converged to a spacetime configuration. Physical consistency and generality imply 
that it is the physically unique one that contains the initial data slice. We also note that, 
as in the situation in Figure |31 the convergence is eventually lost in some simulations; these 
simulations are no longer solving Einstein's equations. 



It is of interest to ask why our approach has not been implemented previously. A number 
of factors were at work. It has been universally assumed that the computational overhead 
of elliptic solvers is excessive. Choptuik indicates a time penalty of x2, incurred by a 
fully constrained 2d evolution, compared to free evolution. This cost is justified by the much 
longer physical lifetimes achieved in the {2d) evolutions of Ref [2^; Choptuik's factor of two 
in time is considered a small penalty. However, in previous implementations of constrained 
{2d) evolution, one additionally had the problem that the solvers were restricted, for instance 
to conformally fiat situations. This required strong gauge constraints on the evolution, and 
meant that generally one re-solved a strongly nonlinear equation (comparable to the initial 
value problem), on each time step. But we have found, even working with straightforward 
package solvers, that the penalty for our approach is only of order 30%. This is because our 
constraint solver (developed for Kerr-Schild superposed initial data), is in fact completely 
general with no form restrictions on the background. (There are, of course, general conditions 
on the elliptic equations to allow their solution ^4\, but we have encountered no difficulties 
in working with physically realistic configurations). 

Thus our elliptic solver can use backgrounds that are strongly nonfiat. They may be 
strongly nonfiat, but they are already very close to constraint solution. This is so because 
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they arise from evolving one timestep with an accurate time integrator; we use a package 
integrator which is typically fourth order accurate in time, and the spatial discretization is 
also fourth order. There is thus only a very small correction arising from the elliptic solve, 
and minimal outer iteration is required. Further, we take note of the fact that our explicit 
time integration has a certain inherent order of truncation error. Thus we do not in fact set 
the residual limit for the elliptic solve anywhere near machine accuracy. Instead we set it so 
that it produces errors which are consistent in size with the evolution truncation error. This 
reduces the internal iteration in the solver to a very small number. Finally, computational 
resources are now becoming adequate for constrained 3-d evolution. ITfiop/sec Computers 
are now accessible, making this work plausible. It will always be the case that constrained 
evolution is more computation- and memory- intensive than unconstrained, but the time has 
arrived that interesting constrained evolutions are possible. Further, ongoing computational 
infrastructure improvements will make this level of computation generally accessible. Well 
within a decade, desktop access to 10 Tflop/sec will mean that constrained evolutions at 
LISA- or LIGO- relevant resolutions will become unexceptional. 



VIII. DISCUSSION 



We have demonstrated that the analytical formulation is not critical to long term evo- 
lution of single black holes, if the correct constraint subtraction is used. Thus, "modern" 
approaches that pose explicitly hyperbolic approaches are not essential; "traditional" g - K 
methods produce comparably long evolutions. The evidence seems to be that there are many 
formulations and subtraction schemes that lead to long-term single black hole stability; we 
have found an especially simple one. 

However, we have found, consistent with theoretical estimates, that the precise subtrac- 
tion is critical (a fraction to two or three decimals). To attack this problem, we have carried 
out periodic solution of the elliptic constraint equations as part of the time integration, to 
enable fully constrained evolutions of the Einstein equations. Our initial results demonstrate 
dramatic improvement of long term stability of a nonspinning black hole simulation. Because 
we solve the constraints, constraint subtraction is irrelevant in this case. We are beginning 
exploration of the constrained evolution approach in spacetimes involving single moving, and 
multiple interacting black holes. We find substantial improvement from constraint solving in 
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every simulation, but we have not achieved infinite-hved g — K simulations, even with den- 
sitized lapse (which is known to improve the hyperbolicity of the system of equations). We 
have begun an approach where the inner excision is made at a constant coordinate surface 
that coincides with the apparent horizon (essentially spherical, or spheroidal coordinates) 
near the excision region j25|. This approach addresses concerns about the validity of the 
stair-step (LEGO) excision region ( 20|). (Additionally the best long-lived isolated black 
hole simulation to date has been carried out by Ref Q with a pseudo-spectral method with 
apparent horizon conforming coordinates.) 

Our approach uses coordinates which are spherical near the hole, with a cartesian region 
further away. It is hoped that this approach may exhibit some of the good properties and 
long lifetime of the codes described in 
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Unconstrained Black Hole using ADM 

FuU 3+1 simulatiocis 
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— No siibtmctiou 5 pts/M 

— Wilh subtiaction 5 pis/M 

— No subtraction 7.5 pis/M 

— With subtraction 7.5 pts/M 
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FIG. 1: The log of the rms norm of the Hamiltonian constraint violations for constraint-subtracted 
and unsubtracted nonspinning black hole simulations with excision. The simulations were per- 
formed at resolutions of M/5 and M /7.5 on a domain size of itlOM. The long-lived runs employed 
optimal constraint subtraction (Eq. (|16j) ). The short-lived run employed no subtraction. 
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Unconstrained Black Hole using ADM 

Full 3+1 simulations 




Time ( M ) 



FIG. 2: The I2 norm of the Hamiltonian constraint violation for constraint-subtracted and unsub- 
tracted nonspinning black hole simulations with excision performed at a resolution of M/5 on a 
domain sizes of itlOM, itlSAf , or ib20M. All cases employed the optimal constraint subtraction 
(Eq. (|16j) ). As the computational domain size increases, the simulation is increasingly shorter lived. 
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FIG. 3: The rms norm of the Hamiltonian and momentum constramts for simulations of a 
Schwarzschild (nonspinning) black hole with excision. None of the simulations used any constraint 
subtraction. The simulations were performed at resolutions of either M/5 or M/7.5 on a domain 
size of ib5M. The long-lived runs employed constrained evolution as described in the text. The 
short-lived runs were unconstrained. 
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FIG. 4: The log of the rms norm for the Hamiltonian constraint in simulations evolving a nonrotat- 
ing Kerr-Scliild black hole with a spatial domain of [— 10A/...10M] at three resolutions: -/Vf/5, M/7.5, 
and M/10. The excision radius was 0.5M in all constrained and unconstrained cases. All simula- 
tions used a densitized lapse with p = ^. The constraints were solved in the constrained evolution 
cases every 0.05M everywhere on the domain except those points where r < 2.0M. Independent 
residual evaluations for the constrained cases with resolution M/7.5 and M/10 are found in Figure 
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FIG. 5: The rms norm of the diagonal spatial components of the Einstein tensor for the M/7.5 
and M/10 constrained simulations presented in Figure |^ 
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FIG. 6: The z = plane of the Einstein tensor component G, 
constrained evolution presented in Figures lU El 
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